clc;
clear all;
% use section4_35;
% exclude a large share of the zeros contributions to speed up the
% bootstrapping procedure
% The 3.5% sample has all nonzero observations but exludes a share of zero contributions such that the sample
% contains 3.5% of obs. 
%%
load('section4_35.mat');
format short g
%% 
% unpack data needed
donation=d.donation; default=d.default;
relfreq=d.relfreq;
%set seed
rng(56789123);
%%
% exclude AD treatment for analysis
donation= donation(default~=0);
default= default(default~=0);

n10= length(donation(default==10));
n20= length(donation(default==20));
n50= length(donation(default==50));

for i=1:1000;
% draw samples from treatments
don10 = randsample(donation(default==10),n10, true);
def10=10*ones(n10,1);
don20= randsample(donation(default==20), n20, true);
def20=20*ones(n20,1);
don50= randsample(donation(default==50), n50, true);
def50=50*ones(n50,1);


don=vertcat(don10, don20, don50);
def= vertcat(def10, def20, def50);

% calculate Delta: gain from opting out
Deltas=Delta_fun(don, def);
%
likelihood=@(lambda)totallikelihood(lambda, don, Deltas, def, relfreq, 1, 1);
opts=optimoptions(@fmincon,'Algorithm','interior-point', 'FunctionTolerance',1e-12, 'StepTolerance',1e-8,'AlwaysHonorConstraints','bounds');
problem = createOptimProblem('fmincon','x0',[0.8; 0.005; 0.005], 'objective', likelihood, ...
  'lb', [0 ; 0; 0], 'ub',[1; Inf; Inf],'options', opts);
gs = GlobalSearch;
[x1] = run(gs,problem);



lambdas(i, 1:3)= x1;
end

save('lambdas_bootstrapping_1000_global.mat', 'lambdas');


